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I. INTRODUCTION 



The Coulomb operator l/ri2 in electronic Hamiltonians is the ultimate source of all 
technical and computational difficulties in quantum chemistry. If it were absent, the many- 
electron Schrodinger equation would immediately separate into one-electron equations and 
its solution would be straightforward, even for large systems. Recognition of this has led to 
much research on the construction and investigation of methods that replace the Coulomb 
operator by approximations in which the two-electron interaction is, in one sense or another, 
"factorized" . 

One of the oldest approaches is to recognize that, when the two particles are far apart, 
their interaction can be modeled accurately by a multipole expansion. Although this ex- 
pansion originated in classical mechanics^"-, it can also be exploited in quantum mechanical 
calculations and this led to the Continuous Fast Multipole Methods^"-. 

In an alternative strategy, the electron density can be modeled in a relatively small 
auxiliary basis, leading to the so-called density fitting techniques'^"— . It turns out that it is 
generally preferable to model the electric field^"— or potential^ of the density, rather than 
the density itself. Cholesky factorization^Si"— is another member of the density-fitting family 
and, by decomposing the two-electron integral matrix, achieves a fit of the electric field of 
the electron density. 

In a series of papers,— "-^P Hackbusch and co-workers have designed schemes for construct- 
ing tensor factorizations of many-electron objects (including the Coulomb operator) and 
such techniques have recently yielded impressive results in Hartree-Fock^ and correlated^^ 
calculations on a variety of small molecules. 

Pursuing a different tack, other researchers have sought to explore the benefits of partition- 
ing the Coulomb operator into a short-range-but-singular part and a long-range-but-smooth 
part. This idea was introduced to theoretical chemistry by Ewald'^^ in 1921 but has enjoyed 
a renaissance during the past 17 years^"— . The Ewald partition 

1 _ erfc(a;ri2) _^ erf(a;ri2) 

has become especially popular within DFT and a host of functionals have been devised^"— 
to treat the short-range {i.e. the complementary error function) part of the operator. 

In work related to that of the Hackbusch group, we have published a series of papers on 
resolutions of two-electron operators into sums of products of one-electron functions. The 
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first four papers^ — focus on the Coulomb operator and the fourth^ uses a Bessel identity^ 
to show that, if ri + r2 < 27r,~ then 

^ oo 

0nim(ri)0n«m(r2) (2a) 

r 12 , 

nlm 



^nlm 



■v) = 2^2 -6nflji{nr)Yi^{v) (2b) 



where 5ij is the Kronecker deha function, ji is a spherical Bessel function and Yj^ is a 
complex spherical harmonic— 

The fifth paper in the series^ focuses on the long-range part of the Ewald partition ([1]) 
and shows that 



oo 



erf(a;ri2) ,* / n , / x /g ^ 
= >^0„im(ri)0nim(r2) (3a) 

'12 ■^"^ 

nlm 

Kimir) = ^\fK^ ji{2rinUjr) yim(r) (3b) 

where the hn and rjn are the weights and (positive) roots of standard Gauss-Hermite 
quadrature^. We have found that this resolution converges well, and is particularly ef- 
fective for the calculation of long-range exchange energies. When used in conjunction 
with a DFT functional for the short-range part of the Ewald partition it offers a potent 
computational tool. 

In the present paper, we will assume that 



■'nlm 



r) = Qn Ji{Kr) RimiKr) (4) 



where J7i(t) = t~^ji{t) and Rimij) = ''"'^m(r) is a solid harmonic.—'^ Equations pbll and 
( l3bl) are special cases of (jlj). For brevity, we suppress the index on henceforth. 
Resolutions ([2]) and ([3]) convert the two-electron integral 

(ab|cd) = JJ x:(ri)x;(ri)T(ri2)Xc(r2)Xd(r2) dn dr2 (5) 

into the sum 

oo 

(ab|cd) = ''^^{ah\nlm){nlm\cd) (6) 

nlm 

of products of one-electron auxiliaries 

{ah\nlm) = j Xa(r)Xb(r)</'nim(r) dr (7) 



In practice, Eq. (|6]) is not used to calculate two-electron integrals explicitly; rather, it is 
substituted into quantum chemical expressions and then facilitates a variety of cost-saving 
factorizations. For example, in three of the earlier papers^i^i^ in this series, we derived 
expressions for the Coulomb and exchange energies in terms of the {ah\nlm) rather than the 
(ab|cd). Equation (E]) has the same formal structure as found in Cholesky decomposition or 
density-fitting schemes, and hence the integral factorization resulting from the resolutions 
offers similar computational benefits. Moreover the resolution approach avoids the cost of 
an (explicit or implicit) Cholesky decomposition or an (explicit or implicit) density-fitting 
matrix inversion. 

However, the usefulness of resolutions such as ([2]) and ([3]) is critically dependent on the 
accuracy of truncated versions of the sum (|6]) and the ease with which the auxiliaries ([7]) can 
be generated. In investigating the first of these issues, we discovered that surprisingly ag- 
gressive truncations of resolutions provide chemically meaningful Coulomb, exchange, MP2, 
full CI, long-range Coulomb and long-range exchange energies^"— In the fourth and 
fifth papers^!^ of this series, we addressed the second issue by presenting explicit formulae 
for the integrals ([7]) where Xa and Xh are s- or p-type Gaussians. However, in order that 
the resolutions be universally applicable, it is essential that an algorithm be available for 
constructing integrals over Gaussians of arbitrary angular momentum. 

The explicit approacb^i^ rapidly becomes complicated and inefficient as the angular mo- 
mentum grows and it is generally agreed that the optimal schemes for forming Gaussian 
integrals are recursive.—"— Since the auxiliaries ([7]) consist of a mixture of Gaussians, spher- 
ical Bessel functions and spherical harmonics, the standard Gaussian recurrence relations 
(RRs) are not immediately applicable. However, we have recently used Ahlrichs' method to 
find an 18-term RR for a general class of Gaussian integrals in phase spaced and we were 
optimistic that an analogous attack would yield an RR for the auxiliaries. 

In Section II of this paper, the sixth in the series, we derive an RR that allows auxiliaries 
to be constructed recursively from simple ingredients. In Section HI, we discuss an algorithm 
for applying the RR with near-optimal efficiency in terms of memory accesses. The efficiency 
and numerical stability of the proposed algorithm are discussed in Section IV. 
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II. RECURRENCE RELATION FOR AUXILIARIES 



If the basis functions in ([7]) are the real Cartesian Gaussians 

[a| = ix- A^r^iy - Ayfyiz - A,)'^^e-°l'-^l' (8) 
[b| = {x- B,f-{y - Byfy{z - B,f^e-^\^-^\' (9) 

where a = (a^;, ay, a^) and b = {b^, by, b^) are vectors of angular momentum quantum 
numbers,— then an RR may be derived using the approach of Ahlrichs^. 
We define Ij = {6jx, Sjy, 6jz) and the scaled derivative 

Id 

Dt = —^ (10) 

^ 2adAj ^ ' 

and the Boys recurrence relation^^ 

[a+l,|=5,[a| + ^[a-l,| (11) 

shows that a higher Gaussian can be formed from two lower ones. By induction, there exists 
an operator 0(a) that generates [a| from an s-type Gaussian,— i.e. 

[a|=6(a)[0| (12) 

and which commutes with Dj. The operator factorizes 

0(a) = d{a^)d{ay)d{a,) (13) 

and it can be shown that 

d{aj)Z = Zd{aj) + a^zd{aj - 1) (14) 

for any function Z that is linear in Aj, i.e. 

DjZ = z and Djz = (15) 

These properties will be important later. 

The (pnim do not depend on A and B, and therefore 

[ab|n/m] = O(a)O(b)[00|n/m] (16) 
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Combining this with ( !TT|) yields 

[(a+ l,)b|n/m] = O(a)O(b)5, [00|n/m] + ^[(a - l,)b|n/m] (17) 

and reveals, as has been noted elsewhere,—"^ that the recursive properties of Gaussian 
integrals depend critically on the first derivative of their fundamental integral.— 
The fundamental integral here is given^^i^ by 

[00|n/m] = Gab Ji{\P) i?z™(AP) (18) 

where 

C = a + /3 (19) 

P = (aA + /3B)/C (20) 

Gab = g„(vr/C)^/^exp(-AV4C - a/3|A - B|VC) (21) 

but, because of the derivative identity 

5, [Ji(AP)] = (-AV2C)P, Ji+i(AP) (22) 

it turns out to be convenient to generalize (|T8|) and consequently (|T71) to 

[00|n/m](^') = (-AV2C)^Gab Ji+p(AP) Pi„(AP) (23) 

[(a+ lj-)b|^H(P) = O(a)O(b)5j'[00|n/m](P) + ^[(a - lj)b|n/m](P) (24) 
Using the Gaussian derivative 

[Gab] = {Pj - A,) Gab (25) 
the derivative (12^ and the solid harmonic derivatives 

dxRlm = C;^ P«-l,m+l ~ C;"^ Rl-l,m-l (26a) 

5j/i?Zm = (c^^ Pz-l,m+l + Pz_i,m_i) (26b) 

dzRlm = Clm Rl-l,m (26c) 

where 



i_ / (/±m-l)(/±m) 2/ + 1 



2/ + 1 

cim = \l{P-m^)^^ (27b) 



one can then show that 



- ^icJ-JOO\n, l-l,m + - c+,[00|n, / - 1, m - 1]^"+^^) 

A 

iS- 

+ -^ic;JOO\n, / - 1, m + + c+,[00|r2, / - 1, m - 1]^^+^^) 

A 

-^Q^[00|n,Z-l,m](^'+^) (28) 

and the [OOlraZmj^^^ are thus closed under differentiation. 

Substituting f l25]) into and applying the commutation relation f lT^ for 0(a) and 
0(b) produces four new integrals of lower angular momentum because the {Pj — Aj) and 
Pj prefactors are linear in Aj and Bj. The commutation relation introduces no additional 
integrals for the remaining terms, as their prefactors do not depend on Aj or Bj. This finally 
yields the 8-term (for j = x ot y) or 7-term (for j = z) RR 

[{si+lj)h\nlm]^P'^ = {Pj - Aj)[ah\nlm]'^P^ + Pj[ah\nlm]^P+'^^ 

+ ^ {[(a - l,)b|n/m](P) + [{si ~ lj)h\nlm]^P+^^} 

+ {[a(b - l,-)|n/m](P) + [a(b - lj)\nlm]^P+^^} 

- %(cr™[ab|n, / - 1, m + 1]^+'^ - clJah\n, / - 1, m - 1](p+i)) 

A 

+ %(cr„,[ab|n, l-l,m+ 1]^^-^'^ + c+^[ab|r2, / - 1, m - 1]^^+'^) 
A 

-^Q^[ab|n,/-l,m](^'+i) (29) 

which allows us to generate [ab | nlm] integrals from the fundamental integrals fl23p . We also 
note that fl29|) cannot be used if A = 0. However, in this special case, the auxiliaries are 
trivial. 

For clarity, we have derived the RR (12^ in terms of complex solid harmonics Rim- How- 
ever, it is straightforward to construct the equivalent RR for real solid harmonics and, for 
optimal computational efficiency, this is preferable to the complex form. Because the struc- 
tures of the real and complex RRs are identical, the discussion in the following sections is 
applicable to both cases. 
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III. ALGORITHM FOR AUXILIARY EVALUATION 



Possessing an RR that relates complicated [ab|n/m] integrals to simpler ones is a good 
start, but it does not constitute a computational pathway until one specifies precisely how 
it is to be used. We now discuss a simple algorithm whose efficiency is close to optimal. 

We define an [ab] class as the set of all integrals that arise when two Cartesian Gaussian 
shells, with angular momenta a and b, are combined with all the (pnim with < n < N, 
< I < L and —l<m< +1. If the superscript p is missing, it implies p = 0. 

The (pnim do not appear explicitly in our class notation and this refiects the fact that all 
steps in the algorithm must be performed for every (pnim- From this point, therefore, we will 
refer only to the bras of the integrals and the reader must remember that our statements 
refer implicitly to the full set of (A^ + 1)(L + 1)^ kets. 

Suppose that we wish to compute an (ab) class, where the round brackets indicate a 
contracted class. If we think retrosynthetically^ and recall that the Head-Gordon and 
Pople (HGP) "horizontal RR" (HRR) 

(a(b + = ((a + l,-)b| + {A, - 5,-)(ab| (30) 

can be applied to contracted integrals,— it becomes clear that it is best to form the [ab) 
class from (eO) classes, where e = a,a + 1, . . . ,a + b. This tactic is beneficial because, when 
constructing classes with 6 = 0, the fifth and sixth terms in the RR (l29l) vanish. 

Next, we note that the RR forms the [(e + 1)0] class from the [eO]^^) and [(e - 1)0] 
classes (where p is unchanged) along with the [eO]'-^'*"^-' and [(e — 1)0] classes (where p is 
incremented). This suggests an algorithm in which the [eO]*^^^ are formed in an outer loop 
over p = a + b, . . . ,0 and an inner loop over e = 0, . . . , a + b — p. The [00]^^^ classes are 
constructed using fl2^ . 

Table [I shows how the algorithm constructs a (dd) class. The [00]*^^^ class in each V 
generation is found via (!23|) . Each of the other V classes is formed, using (!29|) . from classes 
above it in its own generation and the preceding one. The C classes are formed by contracting 
those in the final V generation and, finally, the H classes are formed using fl30l) . 

Having determined the gross structure of the algorithm, i.e. which intermediate classes 
are needed, we next ask how each of these classes can be efficiently constructed. Because 
modern computer architectures are limited primarily by data access, rather than fioating- 
point operations (FLOPs),— we measure the cost of our algorithm in terms of fetch count. 
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TABLE I. Intermediate classes needed to form a (dd) class 



vo 


VI 


V2 


Generation 
V3 


V4 


c 


HI 


H2 


[00] (4) 


[00] (3) 


[00] (2) 


[00] (1) 


[00] 










[10] (3) 


[10](2) 


[io](i) 


[10] 












[20] (2) 


[20] 


[20] 


(20) 












[30] 


[30] 


(30) 


(21) 












[40] 


(40) 


(31) 


(22) 



A. Construction of [00]*^^^ classes 

The [00] classes are defined by ( l23l) but it is important to note tliat tlie {—X'^/2Q and 
Gab do not depend on / or m, the Ji{\P) do not depend on m, and the Rim{^^) depend 
only trivially on n. It is therefore best to pre-compute these quantities and (if and L are 
moderately large) this pre-computation will represent a negligible fraction of the total cost 
of forming an {ah) class. 

The Ji satisfy the backward recurrence relation^i^^ 

Ji{t) = {21 + 3)Ji+i{t) - eji+2{t) (31) 

and the two initial values j7L+a+b(t) and Jh+a+h-At) can be obtained by calling a standard 
Bessel function routine^^. 

The solid harmonics satisfy the forward recurrence relation^*^ 



Rim - y ^^3^ z Ri-i,m - y 2/-3 P-m^ " ^'-^'"^ ^^^^ 
and the reflection formula 

Ri-m = (33) 

and take the boundary values 



Thus, it is efficient to construct the Ru {0 < I < L) using form the Rim {0 < m < I) 
via (132|) and the Ri-m by (|33|) . 



If the factors in (123|) have been thus pre-computed, the assembly of all the necessary 
[00] classes requires only 

Fi = 3{a + b+ 1) (35) 

memory fetches per ket. 



B. Construction of [eOj^^^ classes 

A quick examination reveals that the RR fl29|) is not equally costly for all integrals. 
Specifically, it is cheapest when aj = (for then the third and fourth terms vanish) and 
when j = z (for then the solid harmonic derivative involves only a single term). Since it is 
desirable to minimize the number of integrals on the right-hand side of (12^ that must be 
fetched, we propose the following boolean for choosing the increment direction j when using 
fl29|) to form a member of the [eO]'-^'* class: 



If 


= 1, 


choose j 


= z 


(3 fetches) 


else if By 


= 1, 


choose j 


= y 


(4 fetches) 


else if 


= 1, 


choose j 


= X 


(4 fetches) 


else if Cz 


>2, 


choose j 


= z 


(5 fetches) 


else if By 


>2, 


choose j 


= y 


(6 fetches) 




else. 


choose j 


= X 


(6 fetches) 



Table HT] lists the best j and associated fetch cost for all Cartesian Gaussians [e| with e < 6. 

Because there are (^^^) bras in an [eO] class, the total fetch cost of forming all of these 
classes is 

a+b a+b—p / „\ / i a\ 

p=0 e=0 V / \ / 

where /x, the average fetch cost of (!29|) . is roughly 4. 



C. Construction of (eO) classes 

The reduction of primitive [eO] classes into contracted (eO) classes is achieved by the 
contraction step 

(eO) = EE^-.^'"^^.^^["0] (37) 
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TABLE II. Increment direction (j) and number of fetches (/) for each integral when forming an 
[eO](p) class by ([29]) 



e = 


= 1 or 3 or 


4 




e = 2 or 5 






e = 6 




e 


j 


/ 


e 


j 


/ 


e 


j 


f 


(1,0,0) 


X 


4 


(2,0,0) 


X 


6 


(6,0,0) 


X 


6 


(0, 1,0) 


y 


4 


(1,1,0) 


y 


4 


(5,1,0) 


y 


4 


(0,0,1) 


z 


3 


(1,0, 1) 


z 


3 


(5,0, 1) 


z 


3 








(0, 2, 0) 


y 


6 


(4, 2, 0) 


y 


6 


(3,0,0) 


X 


6 


(0, 1, 1) 


z 


3 


(4, 1, 1) 


z 


3 


(2, 1,0) 


y 


4 


(0,0,2) 


z 


5 


(4,0,2) 


z 


5 


(2,0,1) 


z 


3 








(3,3,0) 


y 


6 


(1,2,0) 


X 


4 


(5,0,0) 


X 


6 


(3,2,1) 


z 


3 


(1,1,1) 


z 


3 


(4,1,0) 


y 


4 


(3,1,2) 


y 


4 


(1,0,2) 


X 


4 


(4,0, 1) 


z 


3 


(3,0,3) 


z 


5 


(0,3,0) 


y 


6 


(3,2,0) 


y 


6 


(2,4,0) 


y 


6 


(0,2,1) 


z 


3 


(3, 1, 1) 


z 


3 


(2,3, 1) 


z 


3 


(0, 1,2) 


y 


4 


(3,0,2) 


z 


5 


(2,2,2) 


z 


5 


(0,0,3) 


z 


5 


(2,3,0) 


y 


6 


(2,1,3) 


y 


4 








(2,2,1) 


z 


3 


(2,0,4) 


z 


5 


(4,0,0) 


X 


6 


(2,1,2) 


y 


4 


(1,5,0) 


X 


4 


(3, 1,0) 


y 


4 


(2,0,3) 


z 


5 


(1,4, 1) 


z 


3 


(3,0,1) 


z 


3 


(1,4,0) 


X 


4 


(1,3,2) 


X 


4 


(2,2,0) 


y 


6 


(1,3, 1) 


z 


3 


(1,2,3) 


X 


4 


(2,1,1) 


z 


3 


(1,2,2) 


X 


4 


(1, 1,4) 


y 


4 


(2,0,2) 


z 


5 


(1,1,3) 


y 


4 


(1,0,5) 


X 


4 


(1,3,0) 


X 


4 


(1,0,4) 


X 


4 


(0,6,0) 


y 


6 


(1,2,1) 


z 


3 


(0, 5, 0) 


y 


6 


(0,5,1) 


z 


3 


(1, 1>2) 


y 


4 


(0,4, 1) 


z 


3 


(0,4,2) 


z 


5 


(1,0,3) 


X 


4 


(0,3,2) 


z 


5 


(0,3,3) 


z 


5 


(0,4,0) 


y 


6 


(0,2,3) 


z 


5 


(0,2,4) 


z 


5 


(0,3,1) 


z 


3 


(0, 1,4) 


y 


4 


(0,1,5) 


y 


4 


(0,2,2) 


z 


5 


(0,0,5) 


z 


5 


(0,0,6) 


z 


5 


(0, 1,3) 


y 


4 














(0,0,4) 


z 


5 


















Totals for e 


= 1, . . . , 6 are 


11,27,42,65,92,124, 


respectively. 






where Ka 


and Kb 


are the dej 


jrees of contraction of 


(a and (b , respectively, and -Da.fca 


and 



-Dfj^fcj are the associated contraction coefficients. 

The total fetch cost of forming all of the (eO) classes from the [eO] classes is 
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TABLE III. Fetch costs to form a (dd) class 



Primitive Contracted 



vo 


VI 


V2 


V3 


V4 


C 


HI 


H2 


3 


3 


3 


3 


3 










11 


11 


11 


11 












27 


27 


27 


6 












42 


42 


10 


36 












65 


15 


60 


72 



D. Construction of the {ab) class 

The target (ab) class can be formed from the (eO) classes using either (130!) or the more 
general TRra relations-^ but, for simplicity, we will consider only the former here. 

Assuming that the internuclear distance Aj — Bj (which is constant for all the kets) is 
pre-loaded, the horizontal RR requires that two integrals be fetched. Using the fact that 
there are ('^^^) ('^^^) bras in an (e/) class, the total fetch cost of forming (ab) from the (eO) 
classes is 

This can be written as an untidy polynomial which is 2nd degree in a and 6th degree in b. 

IV. COMPUTATIONAL COST AND NUMERICAL STABILITY 
A. Computational cost 

The fetch cost for forming an (ab) class is 

Cost = (Fi + F2 + F^)KaK, + F4 (40) 

The first three cost parameters are multiplied by KaKh because they apply to each of the 
primitive bras. 

The results in Table UTt and (135|) . ( 138|) and (139|) can be used to find the cost of constructing 
any [ab) class of interest. To illustrate this, the costs of each of the intermediate classes in 
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TABLE IV. Cost parameters for forming various {ab) classes 



Fetch cost parameters 



Class 


Size 


Fi 


F2 


F3 


Fi 


(ss) 


1 


3 





1 





(ps) 


3 


6 


11 


3 





ids) 


6 


9 


49 


6 





Us) 


10 


12 


129 


10 





(PP) 


9 


9 


49 


9 


18 


(dp) 


18 


12 


129 


16 


36 


ifp) 


30 


15 


274 


25 


60 


(dd) 


36 


15 


274 


31 


168 


m 


60 


18 


511 


46 


270 


iff) 


100 


21 


872 


74 


776 



Table [T] are shown in Table IIIIl and it follows from these that the cost parameters for a (dd) 
class are Fi = 15, F2 = 274, F3 = 31 and F4 = 168. 

The cost parameters for any class up to (//) can be found in the same way, and these are 
given in Table HVl A cursory glance at this Table reveals that the construction of the funda- 
mental integrals (Fi) and the contraction step (F3) are much cheaper than the construction 
of the [eO]*-^'' classes {F2), except when the angular momentum of the bra is low. The cost 
of the (eO) — )■ (ab) transformation {F4) is significant for bras of high angular momentum 
but the fact that this work is performed on contracted, not primitive, bras means that the 
construction of the [eO]*^^-' via (129|) will be the computational bottleneck for most classes 
arising in typical quantum chemistry calculations. 

The structure of our RR fl29|) is similar to that of the vertical RR (VRR) of Head-Gordon 
and Pople^ and, as a result, their fetch and FLOP costs are comparable. Examination 
of Table HV] reveals that the formation of a single, uncontr acted auxiliary integral requires 
approximately 10 fetches, irrespective of its angular momentum. This is comparable to 
the cost of forming a single, uncontracted conventional two-electron integral by the HGP 
algorithirt^ and implies that the use of truncated versions of ([6]) will be competitive with 
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traditional approaches provided that the number of resolution functions (pnim is less than 
the number of basis function pairs XaXh in the system. 

For calculations with heavily contracted basis sets, there exist algorithms^!^ that intro- 
duce contraction early and these offer significantly enhanced efficiency. Although we will 
not explore this extension here, it is straightforward to use the PRISM approach^ to con- 
struct a contracted version of and this offers a cheaper pathway to (ab) classes where 

KaKb > 1. 



B. Numerical stability 

In Section II, we have derived an 8-term RR fl29l) to construct auxiliary integrals of high 
angular momentum. We have advocated that it be used in conjunction with the HRR in 
section III, and we have carefully considered the computational cost of the resulting algo- 
rithm in section IV A. However, we have not addressed the important practical question 
of numerical stability. After all, if our new RR is unstable, the fact that it is cheap will 
be of little real interest. Unfortunately, although the stability characteristics of two-term 
recurrence relations^i^iS^ - including the HRR— i2I - have been studied extensively, the theo- 
retical numerical stabilities of many-term recurrence relations, such as the VRR and our RR 
(1291) . are less well understood. Therefore, having little option, we have resorted to extensive 
numerical experiments. 

We decided to focus on the construction of a [gg] class, using the HRR to transfer angular 
momentum from A to B. We performed a high- dimensional scan, varying the Gaussian 
exponents (a, (3 = 10"^, 10"^ . . . , 10^), the intercenter distance {R = | A-B| = 0, 1, 10, 100), 
and the resolution parameters {n = 0, 1, 10, 100, / = 0, 5, 15, 25, m = 0, g„ = 1), and using 
to generate the 45 [/cs|n/m] normalized integrals, in both double and extended precision. 
Remarkably, the largest absolute error in the double precision estimates of the integrals was 
1.84 X 10-^^ indictating that the RR ([29]) was completely stable for all of the 184 320 
integrals investigated. We conclude that the numerical stability of the new RR probably 
exceeds that of the widely used HRR— i^I and is, in any event, no cause for concern. 
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V. CONCLUDING REMARKS 



In the present article, we have presented a detailed algorithm for the computational 
construction of the auxiliary integrals (ab|nZm) and we have shown that, although such 
integrals involve Gaussians, Bessel functions and spherical harmonics, they can be generated 
recursively at approximately the same cost per integral as traditional two-electron Gaussian 
repulsion integrals. 

This work opens the door to the routine and widespread application of two-electron 
resolutions within Gaussian-based quantum chemistry. We are currently developing a high- 
performance implementation of our algorithm and we will report applications to problems 
of chemical interest in the near future. 
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